  setwd("C:/Users/Jef/Desktop/Applied Statistics")
# setwd("D:/Mijn Documenten/applied statistics - project")
library(qcc)

# opg1 (a)

# --- INDIVIDUAL DATA --- 
# Check whether the process was statistical in control and normally distributed 
vanes 	<- 	read.table("vanes.txt" )
# generate rational subgroups
vanesdata 	<- 	qcc.groups(vanes$V1,vanes$V2) 

#Informal check for outliers, Box and Whisker Plot - individual data
boxplot(vanesdata, horizontal = TRUE, main="Box-and-whisker plot of individual vanes data")

#Check the normality:
#Kernel density plot - individual data:
plot(density(vanesdata),main="Kernel density estimate of individual vanes data",col="red",lwd=3)
#Normal probability plot - individual data:
qqnorm(vanesdata,main="Normal probability plot of individual vanes data",pch=19,cex=1,fg="red")
qqline(vanesdata,lwd=3,col="blue",lty="dashed")

#Standard Normality test - individual data
shapiro.test(vanesdata)

# --- GROUP MEANS ---   
# Matrix to calculate the means and bind them to 'vanesdata' 
means 	= 	matrix(nrow = 20, ncol = 1)
for(i in 1:20) means[i,1] <- mean(vanesdata[i,])
vanesdatam 	<- 	cbind(vanesdata, means)

#Informal check for outliers, Box and Whisker Plot - groupmeans
boxplot(vanesdatam[,6], horizontal = TRUE, main="Box-and-whisker plot of the groupmeans vanes data")

#Check if the group means are normally distributed.
#Kernel density plot  - groupmeans
plot(density(vanesdatam[,6]),main="Kernel density estimate for groupmeans vanes data",col="red",lwd=3)
#Normal probability plot - groupmeans
qqnorm(vanesdatam[,6],main="Normal probability plot of the groupmeans vanes data",pch=19,cex=1,fg="red")
qqline(vanesdatam[,6],lwd=3,col="blue",lty="dashed")

#Standard Normality test - groupmeans
shapiro.test(vanesdatam[,6])

# --- Removed outliers --------------(niet in verslag gezet!) -------
# first we remove group number 8:
vanes2 	<- 	read.table("vanes2.txt" )
boxplot(vanes2, horizontal = TRUE, main="Box-and-whisker plot")
# generate rational subgroups
vanes2data 	<- 	qcc.groups(vanes$V1,vanes$V2) 
plot(density(vanes2data),main="Kernel density estimate",col="red",lwd=3)
vanesqcc 	<- 	qcc(vanes2data,type="xbar")

# now we remove groups 6, 11 and 19:
vanes3 	<- 	read.table("vanes3.txt" )
boxplot(vanes3, horizontal = TRUE, main="Box-and-whisker plot")
# generate rational subgroups
vanes3data 	<- 	qcc.groups(vanes$V1,vanes$V2) 
plot(density(vanes3data),main="Kernel density estimate",col="red",lwd=3)
vanesqcc 	<- 	qcc(vanes3data,type="xbar")
# --------------------------------------------------------------------------

# opg1 (b)

# Capability Analysis: 
vanesqcc 	<- 	qcc(vanesdata,type="xbar")
process.capability.sixpack(vanesqcc,spec.limits=c(0.5005,0.5055),target=0.5030,nsigmas=3)
process.capability(vanesqcc,spec.limits=c(0.5005,0.5055),target=0.5030,nsigmas=3)

LSL		<- 	0.5005
USL		<-	0.5055
1/2*(USL+LSL)

Cp		<-	3.199 
Cpk		<-	2.774 
2*pnorm(-3*Cp)
pnorm(-3*(2*Cp-Cpk))+pnorm(-3*Cpk)

#c)

#e)
#defining the variables
mean 		<-	0.503332 
StdDev	<-	0.0002605242
delta		<-	0.0005
LSL		<-	0.5005  
USL		<-	0.5055 
Cporg		<-	3.199 
Cpkorg	<-	2.774 

#suppose a positive shift
mupos		<-	mean+delta
Cppos 	<-	(USL-LSL)/(6*StdDev)
Cpkpos	<-	min((USL-mupos)/(3*StdDev),(mupos-LSL)/(3*StdDev))

#suppose a negative shift
muneg		<-	mean-delta
Cpneg 	<-	(USL-LSL)/(6*StdDev)
Cpkneg	<-	min((USL-muneg)/(3*StdDev),(muneg-LSL)/(3*StdDev))

#change
Cpkpos/Cpkorg
Cpkneg/Cpkorg

#f)
vanes 	<- 	read.table("vanes.txt")
n		<-	100

x1 		<-	min(vanes$V1)
x100 		<- 	max(vanes$V1)

samplenumber<- 	rep(1:100, each=1)
vanesdata 	<- 	qcc.groups(vanes$V1, samplenumber)

stats.xbar(vanesdata)
xbar		=	0.503332
S 		<-	sd.xbar(vanesdata)

d1 		<- 	qnorm((1-(-0.7+0.5*log(n))/n))
d2		<-	qnorm(1-(5/(n*n^(1/2))))

LL		<-	(xbar-x1)/S
UL		<-	(x100-xbar)/S

d1		<=	LL && LL	<=	d2
d1		<=	UL && UL	<=	d2
UL		<=	d2

p		<-	0.002
LCLbias	<- 	xbar - qnorm(1-p/2)*S*(1+((qnorm(1-p/2))^2+3)/(4*n))
UCLbias	<-	xbar + qnorm(1-p/2)*S*(1+((qnorm(1-p/2))^2+3)/(4*n))

epsilon	<-	0.1
alpha		<-	0.2
LCLexc	<- 	xbar - qnorm(1-p/2)*S*(1 + ((qnorm(1-alpha)*(1/2+(qnorm(1-p/2))^2)^(1/2)) /(n^(1/2)))-epsilon/((qnorm(1-p/2))^2))
UCLexc	<- 	xbar + qnorm(1-p/2)*S*(1 + ((qnorm(1-alpha)*(1/2+(qnorm(1-p/2))^2)^(1/2)) /(n^(1/2)))-epsilon/((qnorm(1-p/2))^2))

vanesdata 	<-	qcc.groups(vanes$V1,vanes$V2)
vanesqccbias<-	qcc(vanesdata, type="xbar", st.dev=S, center=xbar, limits=c(LCLbias,UCLbias))
vanesqccexc	<-	qcc(vanesdata, type="xbar", st.dev=S, center=xbar, limits=c(LCLexc,UCLexc))

